! This represents Q_x-Q_y + N_spin = integer resonance
!
! There is a huge tune shift on the spin clearly visible
!
! This is a program which uses a version of the ALS represented by a flat file
! The sextupoles are turned off
! The symmetry is broken by giving one sextupole a quad component
! The entire acceleration is faked a la Lucio (God help us)
!
! Tracking at the bottom confirms the Stora-Froissart result

program ptc_geometry
use madx_ptc_module
use pointer_lattice
implicit none

character*48 :: command_gino
type(probe) xs0,xs1,XST
type(probe_8) xs
type(layout), pointer :: als
real(dp) closed(6),BETA_XX,BETA_XY,theta0,N0(3),DN0_DX(3),norm,dx,X0(6),XN(LNV)
type(real_8) ray8(6),theta
type(damap) id,one_turn,a
type(normalform) nf
INTEGER MF,I,N,k,mft,pos
TYPE(FIBRE),POINTER:: P
TYPE(integration_node),POINTER:: t,t1,t2
type(internal_state) state
type(damapspin)  one_turn_S,a_S,A_f,A_spin,A_l,A_nl,NORMALSPIN,M_DS,R_DS
type(normal_spin) nf_S
TYPE(spinor_8) N_AXIS
TYPE(spinor) N_AXIS_ave
TYPE(res_spinor_8) N_AXIS_res
real(dp) dr,r,ag,ag0,gamma0,AR,AI,bR,bI,cR,cI,at,dre,di,art,aitt,fac
 
real(dp) x(6),circ,xt(lnv),prec ,rate
logical first
TYPE(work) E_0
REAL(DP) a13,a11
type(c_ray) cray,xcn
complex(dp) xc(6) , bet, Dtune(3)

type(c_damap) c_map,c_spin0,a_cs,as,a0,a1,a2,rot,orb,e_ly,e_r
type(c_damap) ray,id_s
type(c_normal_form) c_n
TYPE(c_spinor) C_N0,n_y,n_r
type(c_taylor) phase(2), nu_spin,nu_spin_par
!-----------------------------------

first=.true.
Lmax = 10.d0
use_info = .true.
 
spin_tune_def=-1     ! (1)
 
call ptc_ini_no_append

CALL  READ_AND_APPEND_VIRGIN_general(M_U,"ALS_FLAT_fast.DAT")

m_U%start%name="ALS "
write(6,*) "Making an node layout"
CALL MAKE_NODE_LAYOUT(m_U%start)

als=>m_U%start
p=>als%start


!!!! circ is the circumference of the ring !!!! 
call get_length(als,circ)

!


ag0=  1.158630183290000D-003*3.2207052d+00/3.4010665d+00

dr=(0.223648040802882d0-0.220712613797765d0)/1.d-6

ag0=ag0-(-0.062244172321728d0-(-0.220712613797765d0))/dr

rate= 2*0.220705155181564d0-6.224417232165955d-002

write(6,*) ag0
ag0=ag0+rate/dr
write(6,*) ag0

pause 8
p=>als%start
E_0=p
write(6,*) 1.d0/e_0%gamma0I,p%ag/e_0%gamma0I
gamma0=1.d0/e_0%gamma0I
 
!!!!!  Make the lattice linear by zeroing all sextupoles !!!!!! 
!!!!!  Set all magnets spin to g=ag0 !!!!!! 
p=>als%start
do i=1,als%n
 p%ag = ag0  
 IF(P%MAG%P%NMUL>=3) THEN
  CALL ADD(P,3,0,0.0000D0)
  if(first) then
!   CALL ADD(P,2,0,0.1D0)  
   CALL ADD(P,2,0,0.2D0)
   CALL ADD(P,-3,0,0.2D0)
   first=.false.
  endif
 ENDIF
 p=>p%next
enddo
 write(6,*) p%ag,p%ag*gamma0
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

closed=0.d0
state=default0+nocavity0  

    CALL FIND_ORBIT(ALS,CLOSED,1,STATE,c_1d_5)
 
call kanalnummer(mf,"spin_result_barber_2qy.txt")

state=default0+only_4d0+SPIN0
XS0=CLOSED
    CALL FIND_ORBIT_probe_spin(ALS,XS0,STATE,c_1d_5,FIBRE1=1,theta0=theta0)

WRITE(6,*) XS0%X(1:6)
WRITE(6,*) XS0%S(1)%X(1:3)
WRITE(6,*) THETA0/twopi
N0=XS0%S(1)%X
write(6,*) " n0 "
write(6,*) n0

state=default0+only_4d0+SPIN0    !+modulation0



!!!! xs0%ac%x(1) and xs0%ac%x(2) will be the (q,p) of the modulating clock
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
call init_all(STATE,2,1)

a_spin_scale=1.d0+(1.d0.mono.c_%nv)

call alloc(id_S,c_map,c_spin0,a_cs,as,a0,a1,a2,rot,orb)
call alloc(c_n)
call alloc(c_spin0)
CALL ALLOC(C_N0)
call alloc(e_ly,e_r)
call alloc(n_y)
call alloc(n_r)
!print_spin=.true.


    call alloc(xs)
    call alloc(ray8)
    call alloc(theta)
    call alloc(N_AXIS_res)
    call alloc(one_turn_S,A_S,A_f,A_spin,A_l,A_nl,NORMALSPIN)
    call alloc(nf_S)
    call alloc(phase)
    call alloc(nu_spin,nu_spin_par)
    ray%n=c_%nv
    call alloc(ray)

!!!! Polymorphic probe is created in the usual manner   
   ID_S=1
   
   XS=XS0+ID_S     !
   
!!!! get spin polymorphic probe after one turn   
CALL TRACK_PROBE(ALS,XS,STATE,FIBRE1=1)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 


!!!! Copy probe_8 into a complex damap because we use the complex package for normal forms
c_map=XS
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
 
!!!! one resonance is left in the map 
c_n%nres=1
c_n%m(1,1)=0   !  the 
c_n%m(2,1)=2   !  the  
c_n%ms(1)=1*spin_tune_def    ! spin tune
!    so the resonance is   nu_x- nu_y - nu_spin = integer
call c_normal(c_map,c_n,dospin=my_true)
prec=1.d-10
! printing all the tunes
write(6,*) c_n%tune(1:c_%nd)

write(6,*) " spin ",c_n%spin_tune

write(6,*) " Resonance tune distance ",2*c_n%tune(2)+c_n%ms(1)*c_n%spin_tune

!  Here I put A_t in phasor's basis!!!    
c_n%A_t=c_n%A_t*from_phasor()   

! m = As a n a^-1 As^-1
c_map=c_n%A_t**(-1)*c_n%As**(-1)*c_map*c_n%As*c_n%A_t
 
write(mf,*) " Normalised map "
call  print(c_map,mf,DEPS=prec)


orb=from_phasor()*c_map*to_phasor() 
 
call c_full_canonize(orb,a_cs,as,a0,a1,a2,rot,phase,nu_spin)

write(mf,*) " Tunes and spin tune "
call print(phase(1),mf)
call print(phase(2),mf)
call print(nu_spin ,mf)

call factor_ely_rest(c_map,e_r,e_ly,n_y,n_r) 
  
call c_n0_to_nr(n_r,C_N0)

write(mf,*) " The Lie exponent of the resonance term "

Write(mf,*) " C_n0%v(1)  < =   C_n0%v(1)-i_*C_n0%v(3) ! coefficient of  1/2(L_x + i L_z) "
Write(mf,*) " C_n0%v(2)  < =   C_n0%v(2)"
Write(mf,*) " C_n0%v(3)  < =   C_n0%v(1)+i_*C_n0%v(3) ! coefficient of  1/2(L_x - i L_z) "
write(mf,*) " "
prec=1.d-6
CALL PRINT(C_N0,mf,PREC)

!!! Here starts some evaluation of the various operators

write(mf,*) " "
write(mf,*) " Here I set the initial ray of the tracking"
write(mf,*) " "

!!! These are in initial coordinates
fac=1.0d0
Xc=0.D0
Xc(1)=0.01D0*fac
Xc(3)=0.03D0*fac
XcN%x=0.D0

write(6,*)
cray=xc
write(mf,'(6(1x,g12.5))') xc(1:4)

write(mf,*) " "
write(mf,*) " Here I set the initial ray of the tracking"
write(mf,*) " "


xcn=(c_n%A_T**(-1)).o.cray


!!!!!!!!!!!!!!!!!!!!!!!!!!  computation of d_ag

do i=1,4
ray%v(i)=xcn%x(i)
enddo

ray%v(5)=1.d0.cmono.5

nu_spin_par=nu_spin-(nu_spin.sub.'0')
nu_spin_par=nu_spin_par.o.ray 



do i=1,4
ray%v(i)=1.d0.cmono.i
enddo
ray%v(5)=nu_spin_par
a11=nu_spin_par.sub.'0'

ray=ray**(-1)


nu_spin_par=ray%v(5)

call print(nu_spin_par,6)


do i=1,4
ray%v(i)=0.d0
enddo
ray%v(5)=-a11

nu_spin_par=nu_spin_par.o.ray
call print(nu_spin_par,6)
!a11=nu_spin_par.sub.'0'
a11=ag0*(1.d0+ (nu_spin_par.sub.'0'))

write(6,*) ag0,a11


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
write(mf,*) " "
write(mf,*) " Evaluate the resonance and tune shift"
write(mf,*) " "

bet=c_n0%v(1).o.xcn


dtune(1)=(phase(1).o.xcn)-(phase(1).sub.'0')
dtune(2)=(phase(2).o.xcn)-(phase(2).sub.'0')
dtune(3)=(nu_spin.o.xcn)-(nu_spin.sub.'0')


write(mf,*) " Numerical value of beta ",bet

write(mf,*) " Numerical value of beta*beta*",abs(bet)**2

write(mf,*) " Numerical x    tune shift ",real(dtune(1)) 
write(mf,*) " Numerical y    tune shift ",real(dtune(2)) 
write(mf,*) " Numerical spin tune shift ",real(dtune(3)) 
write(6,*) " Numerical spin tune shift ",real(dtune(3)) 
pause 555
write(mf,*) " Spin resonance shift ",real(dtune(1)) - real(dtune(2)) + c_n%ms(1)*real(dtune(3)) 
write(6,*) real(dtune(1)) - real(dtune(2)) + c_n%ms(1)*real(dtune(3)) 

write(mf,*) " Floquet coordinates "
write(mf,*) xcn%x(1:2),  xcn%x(1)*xcn%x(2)
write(mf,*) xcn%x(3:4),  xcn%x(3)*xcn%x(4)
write(6,*) " map computed at this position "
write(6,*) " Floquet coordinates "
write(6,*) xcn%x(1:2),  xcn%x(1)*xcn%x(2)
write(6,*) xcn%x(3:4),  xcn%x(3)*xcn%x(4)

n=200000
r=0.02d0
dr=2*r/n   ! rate of change of the approach to the resonance

rate=abs(bet)**2/dr/4.d0

write(6,*)"formula with complex DA", 2*exp(-rate)-1
write(mf,*)"formula with complex DA", 2*exp(-rate)-1

close(mf)
 
call kanalnummer(mf,"spin_result_barber_2qy_new.txt")

call kanalnummer(mft,"spin_barber_2qy.dat")
  

!!!!  Here checking with Tracking


ag0=(als%start%ag*gamma0-r)/gamma0

 !ag=(ag0*gamma0+n/4*dr)/gamma0
a13=ag0
p=>als%start
do k=1,als%n
  p%ag=ag0
  p=>p%next
 enddo
XS0=CLOSED
    CALL FIND_ORBIT_probe_spin(ALS,XS0,STATE,c_1d_5,FIBRE1=1,theta0=theta0)



X0=0.D0
X0(1)=0.01D0*fac
X0(3)=0.03D0*fac

    
XS0%X(1)=XS0%X(1)+X0(1)
XS0%X(3)=XS0%X(3)+X0(3)
XS0%s(1)%x(1)=0.d0
XS0%s(1)%x(2)=1.d0
XS0%s(1)%x(1)=0.d0



a13=0


do i=1,n   !+5000
 
 if(i<=n) then
 ag=(ag0*gamma0+i*dr)/gamma0
   p=>als%start
 do k=1,als%n
  p%ag=ag
  p=>p%next
 enddo
 endif
 
if(a13-a11<0.and.ag-a11>0) then
 write(6,*) a13,ag,a11

!!!! !!!!!!!!!!!!!!!!!!!!!!!
xs1=CLOSED
    CALL FIND_ORBIT_probe_spin(ALS,xs1,STATE,c_1d_5,FIBRE1=1,theta0=theta0)

   ID_S=1
   XS=xs1+ID_S     
CALL TRACK_PROBE(ALS,XS,STATE,FIBRE1=1)
c_map=XS
call kill(c_n)
call alloc(c_n)
!!!! one resoanance is left in the map 
c_n%nres=1
c_n%m(1,1)=0   !  the 
c_n%m(2,1)=2   !  the  
c_n%ms(1)=1*spin_tune_def    ! spin tune
!    so the resonance is   n_y + N_spin = integer
call c_normal(c_map,c_n,dospin=my_true)
! printing all the tunes
write(6,*) c_n%tune(1:c_%nd)
write(6,*) (c_n%tune(1)) - (c_n%tune(2)) + spin_tune_def*c_n%spin_tune 

!  Here I put A_t in phasor's basis!!!    
c_n%A_t=c_n%A_t*from_phasor()   

! m = As a n a^-1 As^-1
c_map=c_n%A_t**(-1)*c_n%As**(-1)*c_map*c_n%As*c_n%A_t
 
write(mf,*) " Normalised map "
call  print(c_map,mf,DEPS=prec)

orb=from_phasor()*c_map*to_phasor() 
 
call c_full_canonize(orb,a_cs,as,a0,a1,a2,rot,phase,nu_spin)

write(mf,*) " Tunes and spin tune "
call print(phase(1),mf)
call print(phase(2),mf)
call print(nu_spin ,mf)

call factor_ely_rest(c_map,e_r,e_ly,n_y,n_r) 
  
call c_n0_to_nr(n_r,C_N0)

write(mf,*) " The Lie exponent of the resonance term "

Write(mf,*) " C_n0%v(1)  < =   C_n0%v(1)-i_*C_n0%v(3) ! coefficient of  1/2(L_x + i L_z) "
Write(mf,*) " C_n0%v(2)  < =   C_n0%v(2)"
Write(mf,*) " C_n0%v(3)  < =   C_n0%v(1)+i_*C_n0%v(3) ! coefficient of  1/2(L_x - i L_z) "
write(mf,*) " "
prec=1.d-6
CALL PRINT(C_N0,mf,PREC)

!!! These are in initial coordinates

Xc=XS0%x

cray=xc
write(mf,'(6(1x,g12.5))') xc(1:4)

xcn=(c_n%A_T**(-1)).o.cray

write(mf,*) " "
write(mf,*) " Evaluate the resonance and tune shift"
write(mf,*) " "

bet=c_n0%v(1).o.xcn

rate=abs(bet)**2/dr/4.d0

write(6,*)"formula with complex DA", 2*exp(-rate)-1
write(mf,*)"formula with complex DA", 2*exp(-rate)-1

 
dtune(1)=(phase(1).o.xcn)-(phase(1).sub.'0')
dtune(2)=(phase(2).o.xcn)-(phase(2).sub.'0')
dtune(3)=(nu_spin.o.xcn)-(nu_spin.sub.'0')


write(mf,*) " Numerical value of beta ",bet

write(mf,*) " Numerical value of beta*beta*",abs(bet)**2

write(mf,*) " Numerical x    tune shift ",real(dtune(1)) 
write(mf,*) " Numerical y    tune shift ",real(dtune(2)) 
write(mf,*) " Numerical spin tune shift ",real(dtune(3)) 
 
 
write(mf,*) " Spin resonance shift ",real(dtune(1)) - real(dtune(2)) + c_n%ms(1)*real(dtune(3)) 
write(mf,*) " map computed at this position "
write(mf,'(1x,i6,2(1x,E15.7))') i,(i-n/2)*dr,xs0%s(1)%x(2)
write(mf,*) " Floquet coordinates "
write(mf,*) xcn%x(1:2),  xcn%x(1)*xcn%x(2)
write(mf,*) xcn%x(3:4),  xcn%x(3)*xcn%x(4)
write(6,*) " map computed at this position "
write(6,'(1x,i6,2(1x,E15.7))') i,(i-n/2)*dr,xs0%s(1)%x(2)
write(6,*) " Floquet coordinates "
write(6,*) xcn%x(1:2),  xcn%x(1)*xcn%x(2)
write(6,*) xcn%x(3:4),  xcn%x(3)*xcn%x(4)

!!!!!!!!!!!!!!!!

endif

a13=ag

 CALL TRACK_PROBE(ALS,XS0,STATE,FIBRE1=1)

 IF(MOD(I,10000)==0) then
  write(6,*) i*100.d0/n ," % done"
  write(6,*) ag,a11
 endif
 IF(MOD(I,100)==1) then
  write(mft,'(1x,i6,2(1x,E15.7))') i,(i-n/2)*dr,xs0%s(1)%x(2)

 endif
  IF(.NOT.CHECK_STABLE) then 
   pause 999
   stop
  endif
enddo

close(mft)
close(mf)






    
    call KILL(ray8)
    call KILL(one_turn_S,A_S,A_f,A_spin,A_l,A_nl,NORMALSPIN)
    call KILL(nf_S)
    call KILL(xs)
    call KILL(theta)
    call KILL(N_AXIS_res)

999 command_gino="OPENGINO"    
  call context(command_gino)   ! context makes them capital
  call call_gino(command_gino)  
 111 command_gino="MINI" 
  call context(command_gino)   ! context makes them capital
  call call_gino(command_gino)  
 



!!!!!!!!! vaguelynecessary baloney
command_gino="CLOSEGINO"
  call call_gino(command_gino)  


 call ptc_end


end program ptc_geometry




!A_S=NF_S%A_T

!!!! factor the "Courant-Snyder" transformation
!call factor(A_S,A_f,A_spin,A_l,A_nl)  ! -->>>   (A_S%m,A_S%s) = (a_f%m, I ) o (I ,A_spin%s) o (a_l%m,I) o (a_nl%m,I)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
